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verse is studied using a previously introduced Lagrangian scheme. An 
' approximate discrete dynamical system is derived, which describes 

the mass agglomeration process qualitatively. Quantitative predic- 
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of the agglomeration, d ~ (t/to)^, is an increasing function of cos- 
mological time t. The major part of the mass hence always lies in a 
region with decreasing mass density. This shows that one-dimensional 
self-gravitating motion is not sufficient to effectively drive structure 
formation in an Einstein-de Sitter Universe. These results are com- 
pared with analogous investigations for the adhesion model (Burgers 
equation with positive viscosity), where the agglomeration is faster, 
and one-dimensional dynamics is effective. We further study the mu- 
tual motion of two mass agglomerations, and show that they oscillate 
around each other for long times, like two "heavy particles". Individ- 
ual particles in the two agglomerations do not mix effectively on the 
time scale of the interagglomeration motion. 

PACS numbers: 02.50.Ey, 05.60.-k, 

Keywords: Self-gravitating dynamics; Mass density; Adhesion model; 
Burgers equation 

1 Introduction 

Structure formation in the Universe is a rich and challenging problem, which 
touches many sides of Physics. It is generally believed that the presently ob- 
served large scale structures have been generated by the process of the gravi- 
tational instability, acting on initially small density perturbations. The tem- 
perature fluctuations on the uniform Cosmic Microwave Background (CMB) 
radiation provide an image of these fluctuations, which were the seeds of 
the present large scale distribution of matter [6, 7, 2]. The observed CMB 
angular power spectrum is dominated by a peak at 1 degree of arc, and also 
shows structures at smaller scales. 

According to present consensus, the constituents of our Universe are or- 
dinary matter (5%), dark matter (25%) and dark energy (70%). The three 
forms can be summarily classifled by how their energy densities change with 
a cosmic scale factor a: ordinary and dark matter behave as a~^, radiation, 
which interacts with ordinary matter, behaves as a~^, while dark energy, 
at least in the simplest models thereof, is independent of a. Hence, even if 
dark energy dominates today, dark matter and ordinary matter dominated 
earlier. Further, the Universe is flat, in agreement with an early epoch of 



2 



inflation. Within that scenario, the primordial density perturbations were 
random Gaussian with power spectrum P{k) ~ Ak , which is also consistent 
with the data. The ordinary matter content can be deduced from observa- 
tions at e.g. visible or radio frequencies. The presence of dark matter can be 
inferred from the observed dynamics of cosmic objects, particularly from fast 
rotation of hydrogen clouds far outside the luminous disc of spiral galaxies, 
as well as high- velocity dispersion of galaxies in clusters [29]. Dark energy, 
or quintessence [5], is equivalent to a non-zero cosmological constant. A, in 
Einstein's equations [29], and there is recent support for a non-zero A also 
from redshift observations. All of these direct measurements can be com- 
pared with theoretical cosmology and the observed angular structure of the 
CMB. A sequence of peaks should indeed arise from coherent acoustic oscil- 
lations in the baryon-photon fluid during an early epoch. Their amplitudes 
and relative positions provide another series of tests of cosmological models, 
and put a different series of constraints on the parameters of such models, 
c.f. [26] for a recent review. 

The growth of small initial perturbations to large-scale structures is essen- 
tially a problem in classical gravitational physics, where all the above mainly 
enters in the initial conditions, and the how the Universe expands as a whole. 
The problem is nevertheless far from trivial. A principal difficulty is that 
a self-gravitational medium has no ground-state, and the dynamics is not 
ergodic. Contrary to e.g. molecular dynamics, one cannot appeal to a shad- 
owing lemma, and argue that even if a simulation only solves the equations 
of motion in a rough approximation, the simulated system nevertheless stays 
close to the real system, with different initial conditions. The single most in- 
fluential result in structure formation, Zeldovich' pancake theory, was indeed 
developed without any recourse to simulation at all, but purely by extending 
linear analysis into the non-hnear domain. The fact that this model and the 
adhesion model are still the benchmarks today shows the difficulty of firmly 
establishing the detailed characteristics of the mass agglomeration process. 

The goal of this paper is to investigate the non-linear regime of the gravita- 
tional instability in the special case that the initial perturbations are planar. 

We also assume an Einstein-de Sitter Universe, and hence disregard the dark 
energy component. With these limitations, one can introduce a Lagrangian 
integration scheme, called the Quintic (Q) model [1, 10], which is both fast. 
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and exact from collision to collision. It is worth stressing that the Q model 
is just one of the possible representation of the dynamics [10]. A concep- 
tually equivalent formulation was previously proposed by [27]. These two 
approaches share the characteristic that they allow to avoid sources of addi- 
tional errors, like truncation, involved in any approximate numerical scheme 



The focus is here on the inner structure of a pancake. The novelty with 
respect previous investigations [8, 27] is that, by making use of a discrete map 
to approximate the exact Q dynamics, we derive a close analytical expression 
for the density profile. These results are validated by direct numerical sim- 
ulation, and compared with analogous results for the adhesion model. We 
further investigate the pair-wise motion of two pancakes under their mutual 
attraction. 

The paper is organized as follows: in Section 2 we present the general back- 
ground, and in Section 3 we re-derive, for completeness, the Q model. In 
Section 4 the evolution of a two particle system is considered: an approx- 
imate discrete system is derived and validated by numerical investigations. 
The analysis is extended in Section 5 where the inner structure of a single 
cluster is investigated. Improving previous results in [10] we establish the 
density profile of a collapsing cluster. In Section 6 we study the dynamics of 
two massive clusters, and in Section 7 we sum up and discuss our results. 

2 From Vlasov-Poisson equations to the ad- 
hesion model 

The evolution of collision-less dark matter in a three dimensional expanding 
Universe is described by the kinetic Vlasov-Poisson equations. Assume an 
inertial reference frame, and label the position with r. Then it is customary 
to introduce the co- moving coordinate x by the following transformation [21]: 



where the scale factor a{t) is function of proper world time. For an Einstein- 
de Sitter Universe 



(see [8]). 



r = a{t)x , 



(1) 




(2) 



4 



where to'^ = QnGpitQ), pito) is the homogeneous density at time to [21, 29] 
and G is the gravitational constant. The Vlasov-Poisson equations then read: 

dtf + ^- Vx/ - V^V' • Vp/ = 

(3) 

VV = 47rG'a2(p-p,) , 

where p is the variable conjugated to x; /(x, p, t) is the distribution function 
in the six-dimensional phase space (x, p) ; is the gravitational potential 
and pfy is the mean mass density. The particle density p(x, t) and velocities 
u(x, t) are given, in term of /(x, p, t), as: 

TJl f 

P{^,t) = — f{^,P,t)dp , (4) 
d J 

p(x,t)u(x,t) = j Pf{^,P,t)dp. (5) 
It is well known, see e.g. [28], that (3) admit special solutions of the form 

/(x, p, t) = ^-P^S'ip - mau(x, t)) , (6) 
m 

where d is the dimension of space and 6'^{.) the d-dimensional delta function. 
We will refer to this class as to single-speed solutions, because to each given 
(x, t) corresponds a well defined velocity u. Assuming (6), a closed system 
on the hydrodynamic level can be derived using (4) and (5): 

dtp + 3-p + -V • (pu) = 
a a 

atu + -u+-(u- V)u = g ('^^ 
a a 

_ V • g = -47rGa{p - pb) , 

where we have introduced g = —Vi/j/a such that V x g = 0. It should be 
stressed that system (7) is vahd as long as the distribution function /(x, p, t) 
is in the form (6). Beyond the time of caustic formation, when fast particles 
cross slow ones, the solution becomes multi-stream. Hence, the pressure-less 



5 



and dissipation-less hydrodynamical equations of (7) are incomplete. 

An ansatz that permits further progress is the so-called condition of par- 
allelism which requires that the peculiar velocity is a potential field, that 
remains parallel to the gravitational peculiar acceleration field [21, 28, 4]: 

g = F{t)u (8) 

The important assumption, in the present discussion, is that the proportion- 
ality F{t) is space-independent. That can only strictly be true in the linear 
regime, where 

F{t) = AnGpbb/b , (9) 

and b is the amplitude of the perturbation. A natural further assumption is 
that only the linearly growing mode is excited [21, 4], which gives 6 ~ in 
our case. Defining the new velocity field v = u/(a6), system (7) reduces to 
free motion in Eulerian coordinates 

dbv + (v V)v = , (10) 

This is Zeldovich' pancake model in the original formulation [25]. It is also 
to referred as the Zeldovich approximation, because the contribution of the 
decaying mode of the density field has been neglected. The ansatz that the 
condition of parallelism holds also in the multi-stream region could be ap- 
proximatively true, but is still a bold assumption to make. 

After caustic formation one may think that the resulting change in the gravi- 
tational force could be modeled by an effective diffusive term. In the adhesion 
model, [15, 13], one hence introduces a term of the form i/V^v in the right 
hand side of the equation (10): 

' dbV + (v • V)v = i/VV 

< v = -VVi (11) 

where ip = il)/{bF{t)). In order for the diffusion term to have a smoothing 
effect only in those regions where the particles crossing takes place, the phe- 
nomenological viscosity parameter u should be small. 
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Although numerical experiments suggest qualitative agreement, the exact 
relationship between (3) and (11) is still an open problem. In a recent paper 
[10] a derivation of (11) was outlined in one dimension. Comparisons with 
numerical results suggested however that at least the diffusion term should 
be modified, to give instead transport equation in the class recently studied 
by Buchert and co-workers [3, 4]. We will here study the system (3) directly, 
and establish the inner structure of a mass agglomeration without the ansatz 
of parallelism. 

3 The Q model, and its relation to the Zel- 
dovich approximation and the adhesion model 

We will consider the evolution of the one- dimensional perturbation in the ex- 
panding three dimension Universe, using the discrete approximation of initial 
condition. Let us assume, that in some interval L along the axis x all the 
matter is concentrated in the N sheets, which we will from now on call par- 
ticles. 



The Newtonian equations of motion for such N particles, interacting via 
gravity, follow from the Lagrangian [21]: 

^ = Y^ - m(f>iri, t) , (12) 

where = Ai:Gp. In the point particle picture the density profile reads: 

p{xi,t) = ^mja~^6{xi - xj) , (13) 

Xj 

where Xi is the co-moving coordinate of the i'th particle, in the direction of 
which the density and velocities vary. 



Expressing (12) as a function of the proper coordinate, Xi, and assuming 
(13), the equation of motion of the i'th particle reads: 



d Xi a, dxi 



+ - 'i'KGpb{t)Xi = a-^Egrav{Xi, t) , (14) 



dt"^ a dt 
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where pb{t) is the mean mass density at time t and 

Egrav{xi,t) = -2'KG'^mjSign{xi - Xj) . 



(15) 



In an Einstein-de Sitter Universe we can make a nonhnear change of variable, 
T = tQlogt/to, such that (14) takes the form: 



dF'x. 



1 dxi 



-\ i Xi — E, 

dr^ Sto dr 3tl * 



grav 



y^ii 'T) 



Q model , 



(16) 



where t^"^ = 67rG'p(to)- This is the representation of the Vlasov-Poisson equa- 
tions which we call the Quintic (Q) model. 

The interest of this formulation is that, as for the classical static self-gravitating 
systems in one dimension, Eg^av is a Lagrangian invariant, proportional to 
the net mass difference to the right and to the left of a given particle, at a 
given time. Hence, in between collisions, (16) has an explicit solution [1]: 

x,{t) = c\ exp( \^ ^ ) + 4 exp(-^^-^) + , (17) 



3^0 



to 



where — —{3tl/2)Egrav{xi,T), is constant. The coefficients c{ and c\ are 
determined by = a;j(r") and = a;j(r"'), i.e. by the states of the particle 
at the time of the last crossing, and read: 



- [x^ + to< - K,] 



(18) 



X, 



-towl - K, 



The form of equation (17) suggests introducing an auxiliary variable z — 
exp((r — r")/3to)- The crossing times between neighboring particles (i.e. 

+ can, hence, be computed by solving numerically the following quintic 
equation: 



in 



0, 



(19) 
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where: 



A. 



Ax^ + toAw^ - {Kl^, - K, 



-(Xf+i - K^) = const 



(20) 



2 

5 

and Ax" = x'^j^^ — x^, Awf = w^^^ — w^. Thus, the evolution of the system 
is recovered by using a version of the event-driven scheme discussed in [20] 
and [11]. Details of the numerical implementation are given in [10]. 



Let as assume that in the interval L all the particles have equal mass rrij 
Lpo/N. Then, the gravitational force in equation (16) transforms in 



2,Ntl 



^sign(a;j - x^) 



L 



2>Nti 



(21) 



where — Ni^right ~ N^^ft is the difference between the number of particles 
to the right and to the left of a given particle. At the initial time r = r° = 0, 
the positions of all the particles are ordered (xj+i > Xj), and we have — 
N — {2i — 1). With in the solution (17) constant we get 



K° = -—N^ = -- + —ii - -) 
2N ' 2^N^' 2> 



(22) 



If we now pick the specific initial conditions such that all the c^^s in (17) 
are zero, we recover, in the discrete setting, the special case of correlated 
velocity and density perturbations, such that only the linearly increasing 
mode is excited. Note that this separation is valid until the first particle 
crossing, i.e. well into the non-linear regime. In this time interval we have 



x.(r) = -^,toexp(-) -- + -(.--), 



(23) 



Hence, defining the new "time" h = 6o exp(2r/3to) (with density dimension), 
the solution of (23) reduces to free motion of the particles. In other words, we 
have re-derived Zeldovich' approximation in the framework of the Q model. 
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In the Zcldovich approximation, this solution is assumed to hold for all time, 
even when the solution becomes multi-stream, while in the adhesion model 
we have merging of the particles after their crossing. The real system must 
necessarily lie somewhere in between. 



4 A two particles system: an iterative map 

Consider two particles of equal mass m interacting by the Q dynamics, con- 
fined in a box of size L. We choose units of length and mass such that 
L — 1 and m = 1/2. To streamhne the following, we first first introduce the 
dimensionless variables: 

^ = ^, (24) 
Q^{0) - J, (25) 
m - (26) 

Thus, before the time of first crossing, equation (17) takes the form: 

<2f) 1 
Qiie) = c{ exp(-) + 4 exp(-^) ^ - . (27) 

where c\ and are given by (18). We select a special class of initial condition: 
the velocities of the particles are arbitrarily assigned, while the positions are 
determined by setting C2 to zero (23). This simply means that there is no 
decaying mode initially, as discussed in the end of the previous section, and 
we can therefore immediately integrate the equations of motion of every par- 
ticle up to its first crossing. In addition we chose a system of reference such 
that the center of mass is at rest at the origin. That means that the initial 
velocities are = —W2 > 0. For simplicity we will suppress in the follow- 
ing the index i and consider the evolution of the leftmost of the two particles. 

The time of first crossing, 9cross, is then deduced from equation (27) and 
reads 

9,ross^ -h^PiO). (28) 
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The rcscalcd particle's velocity is (3cross,o = /5(^cros.s) = 1- Remark that 
(Across is solely determined by the initial particle velocity, while Pcross,Q is a 
dimension-less constant, not dependent on the states of the particles at ^ = 0. 
This means that while the time at which the particles cross {0 cross) depends 
on the initial velocity, their further evolution does not. On the other hand, 
after crossing the initial balance between positions and velocities no longer 
holds, and the condition C2 = is no longer true. From that time on both 
the growing and decaying modes in (27) have to be considered in the analysis. 

For times larger than 6 cross, the separation between the particles increases. 
Conversely, as an effect of the gravitational attraction, their velocities are 
progressively reduced. Hence, a turning point is reached where the elon- 
gation is maximal and the velocities are zero. The inversion time, 9turn, 
follows by differentiating equation (17) with respect to time, and imposing 
the condition qiOtum) — 0- Thus: 

^turn ^cross ~l~ ^ (^^) 

where 6 ~ fin 6. The corresponding maximum elongation, qmax,o — Qidturn), 
reads: 

qmax,o = ^ - - — et = a ~ 0.08. (30) 
It does not depend on the initial velocity. 

Note that the particle inversion occurs before the initial position (i.e. \qi{0) \ = 
1(1-/5(0)) is recovered: this is the effect of the friction-hke term in equation 
(16), which is responsible for the progressive localization of the system around 
the center of mass. An analogous damping is shown to occur for the particles 
velocities. In fact, for 9 > 6 turn, the particles move inward approaching the 
time of successive crossing. Due to the inner symmetry, the latter will take 
place in the origin. The particle rescaled velocity corresponding to the first 
return reads: 

l^ret — Pcross,! — —^-SSPcrossfi — —0.83. (31) 

By a simple generalization of the previous discussion, it can be shown that, for 
larger times, the particle executes damped oscillations. The particles move 
back and forth, displaying a progressive reduction of their mutual distance: 
the crossings become more frequent as time is increased. Take an index 
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n to label particle crossings in the system. The variable 6n represents the 
time elapsed between the n — I'th and n'th crossings, and Pcross,n is the 
dimensionless velocity of the particle at the n'th crossing. Then we make the 
further assumptions of smallness: 

|/3cross,ra| 1- (32) 

For large enough n these are reasonable approximation. Further, we assume 
(32) to hold even for an early stage of the evolution: we will return on this 
point, and provide numerical support to the validity of this ansatz. Within 
these approximation the dynamics of (17) is reduced to the following simple 
map: 

Pcross,n+l Pcross,n ^Pcross,nj 

= 2|/?cross,n| (33) 

where /^qmax,n stands for the maximal elongation reached by the particles 
after n encounters. It is interesting that system (33) is actually equivalent 
to the motion of an impurity advected by an unsteady Burgers flow with a 
particular numerical relation between frictional and inertial forces, see [14]. 



Consider then the finite difference formula: 

^l/^crossl |/5cross,n+l| |/3cross,n| 



n+l 



(34) 



and rewrite the right hand-side using the first two relations of (33). Approxi- 
mating the finite differences with differentials and integrating, using as initial 
condition for \l3crossi6 = 9cross) \ the velocity of at first crossing \Pcross,o\ and 
for \Aqmax{0 — 9turn)\ the first maximum elongation (30), one ends up with: 

\Pcross{9) \ = |/3cross,o| ^Xp ^~g(^ ~ ^cross)^ , Pcross,0 = 1; (35) 
IMmaxMl = '^Qmax exp ("^(^ " ^turn)^ , Qmaxfi = Ot. (36) 

Here \Pcross{9)\ is the dimension-less velocity of the particles at the moment 
of their crossing and Aqmax{d) the maximum separation between the two par- 
ticles, at rescaled time 6. Both the velocity and the inter-particle distance 
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(and also the time elapsed between two successive crossing 6n{0) ~ \(^cross{0)\) 
decrease exponentially, as function of time 9, or, equivalently, a power-law 
decays is displayed in term of cosmological time t: \/3cross{t)\ ~ t~^^^ and 

Let us turn to numerical simulations to test the reliability of our results. 
The initial condition is displayed in the small inset of Figure 1: here C2 = 0, 
thus before the time of first crossing, only the contribution of the growing 
mode has to be considered. In the main plot of Figure 1 the absolute value 
of the particle velocity, Across, is represented as function of rescaled time 6, 
in scale lin-log. The solid line refers to the analytical prediction (35). In 
Figure 2 the particle maximum separation Aqmax{()) is plotted versus 9, in 
hnear-logarithmic scale (36). The sohd line shows the theoretical result (36). 
In both case good agreement is displayed, even in the initial stage of the 
dynamics, thus confirming a posteriori the validity of our assumptions. 

5 Density profile: Q model vs. adhesion model 

The objective of this section is to investigate the density profile of a collapsing 
cluster in the Q model, and compare that to the adhesion model. We show 
that collapse in the Q model, which we recall is nothing but one-dimensional 
self-gravitating dynamics on the background of an Einstein-de Sitter Uni- 
verse, is less pronounced. We derive a general expression for the density 
profile. To this end, we begin by analysing the motion a particle in an 
external halo, moving in the potential from a much more massive and very 
localized agglomeration. Then, we reconstruct the whole density distribution 
by combining together the contribution from each particle. These theoretical 
results are then compared first with numerical simulations of the Q model, 
and then with predictions based on the Zeldovich approximation and the 
adhesion model. 

Consider a system of particles of equal mass, confined in box of finite 
size L = 1, symmetric with respect to a reference point located at the origin. 
First, assume N — 21 particles to form an inner bulk of width A, localized 
in the center of the box, with I much smaller than N, and A much smaller 
than L. The remaining 21 particles are assumed uniformly distributed in the 
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outer regions at the borders of the box. 



To first approximation we can neglect the fine structure of the cluster, and 
replace it with an heavy, structure-less, macro-particle, carrying the same 
amount of mass. In addition, the gravitational interaction between the par- 
ticles in the halo is also neglected. As a consequence, these particles arc then 
not affected by the discreteness of the inner distribution: a sudden change 
in the acceleration is experienced when the origin is crossed. 

As in the previous section, we set initial conditions such that only the grow- 
ing mode is active, before the first particles crossing occurs, and use the same 
dimension-less variable. Consider the motion of particle i, situated to the left 
of the agglomeration and define: 

7< = i4(.-i). (37) 

The positive parameter 7j represents the absolute value of the initial position 
of the i'th particle {w^ = in (23)). In analogy with (28)-(31), one finds: 



(^i,cross — —-^^T^ y-^j , (38) 

PijCross Pii^cross) ^li i (^9) 

3 

^i.turn ^i.cross ~l~ ^ ^i,cross ~l~ "r ln 6 , (^0) 

O 

qi,max,o = 4:a-fi . (41) 

In this approximation, the i'th particle crosses for the first time the origin 
with rescaled velocity Pi,cross,o- Note that the latter is solely determined by 
the particle initial position. The corresponding maximum elongation qi^maxfi 
is 4q; ~ 0.32 times smaller then the unperturbed distance 7j. 

By repeating the same procedure as in Section 4, one ends up with the 
following relations for, respectively, \l3i^cross\ and \^qi,max\'- 

\Pi,cross{0)\ = exp (-^(^ - Oi,cross)^ , (42) 
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\^qi,max{0) \ = ia-fi exp 



(43) 



To validate our analysis we performed numerical simulations and studied 
the single particle behavior. In our experiments we considered the following 
initial condition: a finite, but large, fraction of all the particles, say N — 2£, 
is uniformly distributed in a narrow region, say A, centered around the ref- 
erence origin. The remaining particles, I of them, are symmetrically placed 
in the two lateral regions, with spatially uniform distribution. As a conse- 
quence, the inter-particle distance in the central bulk is smaller than in the 
external halos. In the limit when A — > and £ = 1, we reach the conditions 
assumed in the preceding derivation. Velocities are given as a smooth func- 
tion of positions. 

In the upper panel of Figure 3 a phase space portrait of a single particle is 
represented: the damping, both in q and /3, is clearly displayed. In the lower 
panel, the rescaled particle position q is plotted vs. time 6: the thin solid line 
refers to the simulation, while the thick curve is plotted from relation (43). 
Numerical simulations show good agreement with analytical predictions (42) 
and (43). 

The previous results can also be extended to the case when particles are 
initially uniformly distributed in the finite size box. Given particle j, we ne- 
glect the interaction with the external masses (i.e. \qj\ > and mimic the 
attraction of the inner particles by considering the forces exerted by a sin- 
gle massive object, localized in the origin. In this approximation, according 
to the hues of the preceding discussion, the behaviors of the \Pi,cross{^)\ and 
\^Qi,m.axi0)\ are also described by equations (42) and (43). Despite the drastic 
assumptions involved, numerical simulations show relatively good agreement 
with theoretical predictions, provided the velocity distribution are smooth 
enough. 

The goal of the remaining discussion is to derive an analytical estimate 
of the density profile by making use of these results. In the continuum limit 
(N — > oo), the density p(g, 9) reads: 



where we introduced the Jacobian of the transformation from Lagrangian 




(44) 
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to Eulerian coordinates. The density field has a rather complex structure 
including singularities p{q, 9) ~ — 9m)^^^^, where qm stands for the particle 
turning point. A rough characteristics of the cluster size, can be deduced 
from equation (41). In fact, before the time the outmost particle reaches, as 
the last one, the massive bulk, the size of the agglomeration, g^, is given by: 



qM = Aa^,,{e) , (45) 

where — 7i*(^) is the initial coordinate of the particle which falls into the 
origin at the time 9 — 6. From (38), (40) it follows: 

74^) = ^ exp C-{9 - ~9)) = ^ exp ("-9) (46) 

When the continuous limit is recovered, = /?(— 7*), 7* being the La- 
grangian coordinate of the most periferic particle of the cluster, at time 
9. Then assume that P{—^) — > 0, as the edge of the box is approached 
(7 1/2)- Then from (46), 7* tends to 1/2, at large times, and conse- 
quently, from (45), qs{9) 2a; ~ 0.16. This observation on is in good 
agreement with the results of numerical simulations presented in a in previ- 
ous work [10]. 



Inside the pancake a multistream flow has developed: the density at co- 
ordinate q follows by summing the contribution of each stream. Consider the 
mean density distribution assuming the length scale Aq much larger than the 
distance between two successive peaks i.e. larger than \qm+i — qm\- The sum 
can be then approximated by an integral over the Lagrangian coordinate 7. 
By introducing the density distribution function fp{'y, 9) associated to the 
Lagrangian coordinate 7, we can write: 

/"(max 
fph,0)dl. (47) 
,'min 

Assume that all the particles belonging to a small Lagrangian interval, A7, 
are uniformly distributed inside the corresponding Eulerian segment Aqmax{9, 7) 
(43): 

4 

Aw(^,7) = \^q^,maxm = ^xp f--^) , X = A'^^Q'^^'a (48) 
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and consider the mass conservation: 

P0A7 = (/(7, ^) A7) ^qma.{e, 7), (49) 

The interval [7mm(?, ^), 7maa;(6')] selects the particles that contributes to the 
density in g, at time 9. In particular they should have reached the origin, 
thus: 

l^inmax) /2 



exp (3^) = Imax- (50) 



4 

On the other hand the amphtudes of the particles' oscillations, ^qmax{9-,l) 
have to be larger than thus ^min is found by solving the following implicit 
equation. 



/ \ I mm / c 1 \ 



Finally, the density reads: 



PiQ, 0)= Po -r 7T-^ = — exp / —d^, (52) 

Jlmin i^qmaxKyil) X / J Imin 78 

and p(g, ^) = when: 

\(l\>(lw{0) = Mmax{d ) ^max) ■ 

(53) 

Here qw{0) is the width of the mean density distribution. The density profile 
(52) depends on the initial particle velocity, namely /9(7). In the following, we 
will solve analytically the integral in (52) for special class of initial conditions, 
and compare with numerics. In addition, we will discuss the analogous case 
in the adhesion model. 

5.1 Step Profile of Initial Velocities 

Assume the initial velocities profile to be given by: 

P{q) = -/3osign(g), (54) 
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where (3o is a constant coefficient and, without loss of generahty, we focus on 
the region q > 0. 



First consider the evolution of the density in the Zeldovich approximation, 
with initial positions as in equation (23). Since the velocity profile at = 
is a step function, the particles with initial negative (rcsp. positive) spatial 
coordinate will move as a whole towards the origin, keeping their mutual 
distances constant. Thus before the last left particle reaches the origin. 



a two-stream fiow develops over a finite region of width 2\q,, 



w,zeld\ 



where 



I = (3o exp(26'/3)/4. Note that the the density in each flow is equal to 
the initial density po- 



Then follow the development under the Q dynamics later in time. By in- 
serting (54) in (52) and performing the integral, one obtains: 



Op 



1/4 



1/3- 



27r, 



(55) 



where ^max is solution of equation (50) and qp{0) is the spatial scale of inner 
structure 



X 



23/51 



exp 



'-e 



(56) 



Wc recall that qwiO) stands for the width of the pancake, see eq. (53). 
Consider also the mass function M(g, 9): 



(57) 



M{q,e)= / p{U)dt 
Jo 

Inserting eq. (55) and performing the integral one ends up with: 

3/4 o / 1 \ 1/3 



M{q,e)=po 



2 V27„ 



(58) 



In the initial stage of the evolution (3{'jmax) — Po, thus, from equation 
(50): 



Imax = Y ) ■ 



(59) 
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The pancake's density distribution is therefore characterized by two length 
scales: the spatial scale of inner structure, qp{9), which decays exponentially, 
see (56), and the width of the agglomeration which grows exponentially ac- 
cording to: 

= ^ exp ~ 0.1/3o exp . (60) 

It is worth stressing that during the initial stage of the pancake's formation, 
the width qyj increases with time slower than as predicted in the Zeldovich pic- 
ture. From equation (55) it follows that the density in the bulk increases ex- 
ponentially in time-like coordinate 6*, according to p{q, 9) ~ exp(6'/6). 
In other words, the system tends to display a progressively denser core, lo- 
calized near the origin. The mass contained in the cluster is also increasing 
exponentially, mp{6) = 2M{q^,9) = poP exp{29/3)/2, in agreement with the 
predictions of the Zeldovich model. The discussion above applies to the co- 
moving frame, but physical density is measured in the inertial frame; the two 
are related by (1). We will return to this point in Sect. 7. Let us however 
in brief state that the collapse in co-moving coordinates is slower than the 
expansion of the Universe. The physical density therefore on the contrary 
decreases at around almost all mass points, and in an inertial coordinate 
system the agglomerations spread out over time, albeit more slowly than the 
expansion of the Universe whole. 

Equations (58)- (60) suggest that, initially, both the mean density distribu- 
tion and mean mass function are self-similar. In particular: 

^mp{9)M{q/q^{9)). (61) 

In the very late time evolution of an isolated cluster all particles have 
reached the inner bulk (i.e. 9 > — 3/2 ln(/5o/2) ). The box length, L, is 
normalized to one so that, mp{9) — po- Then ^max — 1/2 and equations 
(55), (57) indicate the a self-similar collapse: 

M{q,9) = M{q/qp{9)). (62) 

Now both the density and the mass function are solely characterized by a 
unique length scale, which is both the inner length scale, qp{9), and the outer 



M(g, 9) = 
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length scale, g^. In this case is the coordinate of the out-most particle, 
and the density for q > qp is identically equal to zero. 

Summing up, we have shown that for an initial condition in the form of 
equation (54), the mean mass function M{q,6), displays a self-similar col- 
lapse, in different stages of the evolution. Introducing the rcscaled variable 
X = q/qw{d)-, the function M{q/ qy,{6)) takes the universal form: 



M{x) = 



2x^" - -X 
2 



(63) 



for X = [0, 1], where M(l) is equal to one half. From (63) it follows that the 
matter is mainly concentrated around the pancake center. In particular, one- 
half of the whole mass (i.e. 2M{xo,5) = 1/2) is found to lay in a symmetric 
interval defined by Xq.s = 0.14, while 90% of it, is distributed in a segment 
equal to half the size of entire cluster (xq.q = 0.54). 

To validate our theoretical analysis we now turn to numerical simulations 
and consider the late stage of evolution of an initial perturbation in the 
form (54). The particles are initially uniformly distributed in space. In 
Figure 4 the normalized density profile is plotted at late stage of evolution 
(thin solid line) and superposed to the theoretical prediction (55) (thick solid 
line). The same curves are represented in scale log-log in the left inset: here 
the circles refer to the simulation. No parameter need to be adjusted by 
numerical fitting and the agreement has to be considered satisfying. The 
corresponding phase space portrait is reported in the right inset of Figure 4, 
where the characteristic spiral behavior is clearly displayed. In figure 5 the 
cumulative mass function is plotted as a function of the variable ^ = q/qp. 

In a previous work of two of the authors, [10], we studied the late time 
evolution of an isolated perturbation within the framework of the Q dy- 
namics. In particular we measured the progressive contraction of the inner 
region of the agglomeration, when compared to the overall Universe expan- 
sion. This effect was computed by the width of region, Ag^, that contains 
a fraction /i of the whole mass of the system, centered around the position 
of maximum density. It turned out that the interval Ag^ shrinks in time 
according to a power- law of cosmological time t, with exponent —2/9. An 
heuristic interpretation was also provided in [10]. To make a bridge between 
our present investigations and these earlier findings we define q^ (63) such 
that 2M(g^, 9) = po/i = const . From equation (62) it follows that Aq^ 
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scales proportionally to g^; thus we are led to assume = x^Qp, where is 
a positive coefficient. Equation(56) implies: 



?.-^exp(--^^) (64) 



in qualitative agreement with the results in [10]. As an example, let us focus 
on the case /x = 1/2. In this = 0.14. We performed numerical 

simulations starting with the same class of initial condition as in Figure 4. 
In Figure 6, 51/2 is plotted vs. 9: circles refer to the numerical simulation 
and the sohd hne represents equation (64). The slopes of the two curves 
in figure 6 agree, but the amplitudes do not. Wc can make the difference 
smaller by adjusting some of the hypotheses made in the derivation above. 
In particular, we considered the particles belonging to the cluster to be uni- 
formly distributed in the finite interval /S.qraax given by (48). In reality, the 
particles spend more time on the border. This effect is explicitly modeled by 
assuming: 

%,eff = CQp (65) 

where c is an "ad hoc" numerical factor, larger than one. 

To complete this section, let us compare our theoretical expression (55) 
with the density profile derived in the framework of the adhesion model. 
These results, which we recall in the following, are detailed in [12]. Let as 
consider the stationary solution of Burgers equation : 

v{x, h) = Vst{x) = —U tanh , (66) 

where 5 = f//2z/ is the width of the shock. This solution is also the asymptotic 
solution of Burgers equation for the initial the step-function profile previously 
considered. The trajectories of individual particles x{b) satisfy the following 
equation 

where the field v{x, h) is determined by the solution of the Burgers equa- 
tion. Assuming that v{x, h) is a stationary solution of Burgers equation from 
(66), (67) we have the following expression for the coordinates of the particles 

x{b,y) = 5Arcsinh [sinh(|//5) exp(— ?76/5)] (68) 
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where y stands for the Lagrangian coordinate. From the conservation of 
mass, the foUowing expression for the density is derived: 



po cosh ( f ) 

p{x, h) = ' (69) 

/sinh2(|)+exp (^) 



Solutions (68) do not apply to the intermediate stage of evolution. Thus we 
will focus on the asymptotic behavior, for time b ^ S/U and consider the re- 
gions \x\ > dp, where 6p = 5 exp{—Ub/5). The particles in the agglomeration 
behave according to 

x{b,y)^yeM-Ub/S). (70) 

Thus in the adhesion model a faster collapse is displayed, when compared 
to the exact Q dynamics: in the latter a power-law decay in "Burgers time" 
b, is in fact produced, namely qp{b,y) — y{b/bo)~^^^. As concern the density 
profile, for x < Sp, one gets: 

poexp(^) 

p{x,b) ^ ^LlXi2=, (71) 



l+(f)exp(^^ 



hence maximum value of the density increases exponentially in "Burgers 
time", Pmax(0, &) = poexp{Ub/S). In the interval Sp <^ x <^ 5 the density 
transforms into a time-independent power-law distribution 

p{x,b)=po- = pojj^. (72) 

Thus in the adhesion model we have that the asymptotic density distribution 
for the initial step profile is localized in the regions \x\ < 5 — const, has 
non-intcgrablc time independent power law tails (72) with cutoff scale Sp = 
6 exp{—Ub/6). In the Q model we have instead an integrable power-law 
distribution, see (55). 

When considering a perturbation in a finite box, in co-moving coordi- 
nates, both the Quintic and the adhesion models show the formation of a 
dense structure. However, in the adhesion model, the width qp of the pan- 
cake decreases faster ( Qp ~ exp —b) then in Q-model ( Qp ~ b~^^^ ). 
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These discrepancies suggest that, at least in one dimension, the adhesion 
approach is vahd only as an approximate model of structure formation. This 
observation agrees with the conclusion in [10] where a new transport equation 
was proposed. 

6 Interaction of two massive clusters 

In this section we consider the dynamics of two interacting clusters. As a 
first approximation, one can neglect the inner discrete structure of a cluster, 
by introducing a single macro-particle to mimic its dynamics. Thus, the 
problem is first reduced to the study of the evolution of a two particles system, 
investigated in Section 4. According to this picture, the agglomerations move 
back and forth displaying damped oscillation. This observation agrees with 
previous numerical investigations reported in [9] . Following the discussion of 
the previous section, each massive agglomeration experiences a contraction of 
the inner core: the competition of these two effects determines the late-time 
asymptotics of the system. 

Numerical simulation are performed starting with the initial condition 
displayed in Figure 7 a). The distance between the two centers of mass, 
S, is plotted vs. time 9, in Figure 7 c). The damping tendency is evident 
(thin solid line) and the theoretical prediction (36) (thick sohd hues) agrees 
with the numerical findings. At late time the two clusters still present a 
finite separation, see 7 b). This is in a good agreement with the theoretical 
predictions of Sections 4 and 5.1, both the widths of the agglomerations 
and amplitude of their mutual oscillations decay as the same exponential 
function of time, hence their ratio is a constant number. Therefore, the 
system composed by two interacting clusters, with finite inner structure, 
evolves approximately in a self-similar way. 

7 Concluding remarks 

In this paper we considered the problem of structure formation in a Einstein- 
de Sitter Universe, focusing on planar perturbations. This was done by using 
the Q model, or Quintic model, a Lagrangian representation previously de- 
rived in [1], in which one of the steps is to solve a large number of quintic 
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equations. We underline however that the Q model is nothing but the motion 
of a system of gravitationally interacting particles, on the background of an 
expanding Universe, with convenient choices of time and space coordinates. 
We showed that the Q model can be approximated by an even simpler iter- 
ative map. This observation permitted to establish a connection with [14] 
where the motion of turbulent flows with strong friction was analyzed, and 
to derive theoretical predictions for both the density profile and the mass 
function in a collapsing cluster. These predictions were cross- validated with 
direct numerical simulation. 

The main conclusion of this work is that the inner structure of a Zeldovich 
pancake is different from the adhesion model, with finite viscosity in the 
Burgers equation. That there is a discrepancy is what one would expect, but 
this had nevertheless to our knowledge never been shown before. As a math- 
ematical problem, a qualitative difference is that in co-moving coordinates, 
the size of collapsing cluster goes down as power-law with cosmological time 
in the self-gravitating system, but as a stretched exponential in the adhesion 
model ^. As a physical problem, the difference is even more dramatic: in the 
self-gravitating system the characteristic scale of the density distribution is 
given in co- moving coordinates by equation (56), which decreases, with cos- 

_ 2 

mological time t, as (t/to) . The co- moving system of coordinates is related 

2 

to an inertial system by equation (1), where the proportionality is (t/to)^. 
Combining the two, the characteristic scale of a collapsing cluster hence in- 

4 

creases in time, as d{t) ~ [t/to)^ . Since the density distribution of equation 
(55) is self-similar, the maximum density attained depends on the graininess 
of the initial conditions. We have not investigated the law of the increase of 
the maximum density, but it is at least possible that with continuous initial 
density distributions arbitrarily high densities can be produced, also if mea- 
sured in an inertial system of reference. A more interesting characteristics is 
however the fraction of total mass M where (in an inertial system) the den- 
sity is higher than some cut-off Pcr- That quantity scales as {^ ^ ^(t) ) ■ Since 
d{t) is an increasing function of time in the inertial system, that fraction de- 

_4 

creases as (t/to) ^ . In this sense one-dimensional self-gravitating motion, on 
the background of an expanding Universe, is not effective in bringing about 

^In "Burgers time" the collapse is simply exponential, see subsection 5.1. 
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mass collapse. 



In the adhesion model the characteristic scale decreases, both in a co-moving 
and an inertial system. Here both the maximum real-space density and 
the fraction of mass at a density larger than some cut-off increases. The 
adhesion model is therefore, in contrast, effective in bringing about mass 
collapse in one dimension. In other words, the adhesion model over-predicts 
the clustering tendency, and overestimates the peak density. A Universe in 
which most mass lies in regions where density decreases can never give rise 
to galaxies, galaxy clusters and other significant objects. The analysis given 
here therefore suggests, quantitatively, why a one-dimensional description of 
mass agglomeration is insufficient. 

We have also investigated the mutual motion of two pancakes, and shown 
them to be closely similar to the decaying oscillations of two mutually at- 
tracting "heavy particles" . 
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Figure 1: The absolute value of the velocity Pcross is plotted as function of 
rescaled time 6, in scale lin-log. Circles refer to the simulations, while the 
solid line is the theoretical prediction (35). The initial condition is displayed 
in the small inset. 
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Figure 2: The inter-particle distance \Aqmax\ is represented as function of 
rescaled time 9, in scale lin-log. Triangles refer to the simulation, while the 
solid line is the theoretical prediction (36). 
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Figure 3: Upper panel: phase space portrait of a single particle. Lower panel: 
position q vs. rescaled time 9. The thin solid line refers to the simulation 
while the thick curve represent the theoretical prediction (43). No parameter 
needs to be adjusted by numerical fitting. Here N = 1024, N — 2i = 768, 
A = 0.15. 
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Figure 4: Main plot: normalized density profile p vs. at late stage of evolu- 
tion rescaled position q, starting from a step-profile for velocities. The parti- 
cles are initially uniformly distributed in space. The thin solid line refers to 
the simulation, while the thick solid one represents the theoretical prediction 
(55)('ymax = 1/2). In this case N = 2000, 6 = 23.27, f3o = 4.9 x 10"^ Hence, 
Qp = 0.0183. Left inset: normalized density profile p vs. rescaled position q 
in log-log scale. The circles refer to the simulation. Right inset: phase space 
portrait at 6 = 23.27 . 
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Figure 5: Mass function M{q, 9) in lin-log scale vs. normalized coordinate 
q/qp. The particles are initially distributed uniformly in space, velocity a step- 
function of position. The thick solid line represents the theoretical prediction 
(58)('yrnax = l/2j. Thc circles, stars and crosses refer to the simulation at 
different times. 
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Figure 6: gi/2 is plotted vs. 9. The circles represent the results of the numer- 
ical experiment. The solid line refers to the theoretical prediction ( equation 
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Figure 7: Figure a): Initial phase space portraits. Here N = 200. Figure b): 
Phase space portrait at time 9 = 25.5. Figure c): 6 vs. rescaled time 9. The 
thin solid line refers to the simulation. The thick curves represent eq. (36). 
Here 9rrn-i. = 1.38. 
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